//
// File: pburg.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 20-Mar-2025 09:03:00
//

// Include Files
#include "pburg.h"
#include "arburg.h"
#include "calcAR_data.h"
#include "datawrap.h"
#include "fft.h"
#include "freqz.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include <cmath>

// Function Declarations
static double rt_hypotd_snf(double u0, double u1);

// Function Definitions
//
// Arguments    : double u0
//                double u1
// Return Type  : double
//
static double rt_hypotd_snf(double u0, double u1)
{
  double a;
  double y;
  a = std::abs(u0);
  y = std::abs(u1);
  if (a < y) {
    a /= y;
    y *= std::sqrt(a * a + 1.0);
  } else if (a > y) {
    y /= a;
    y = a * std::sqrt(y * y + 1.0);
  } else if (!std::isnan(y)) {
    y = a * 1.4142135623730951;
  }
  return y;
}

//
// Arguments    : const ::coder::array<double, 1U> &x
//                double p
//                double varargin_1
//                double varargin_2
//                ::coder::array<double, 1U> &varargout_1
//                ::coder::array<double, 1U> &varargout_2
// Return Type  : void
//
namespace coder {
void pburg(const ::coder::array<double, 1U> &x, double p, double varargin_1,
           double varargin_2, ::coder::array<double, 1U> &varargout_1,
           ::coder::array<double, 1U> &varargout_2)
{
  array<creal_T, 2U> den;
  array<creal_T, 2U> num;
  array<creal_T, 1U> h;
  array<double, 2U> a;
  array<double, 2U> b;
  array<double, 2U> b_a;
  array<double, 1U> Sxx_unscaled;
  array<double, 1U> a__1;
  array<double, 1U> wf;
  double d2_tmp;
  double d2scaled;
  double delta1;
  double v;
  int i;
  int k;
  int nx;
  arburg(x, p, a, &v);
  d2_tmp = varargin_2 - varargin_2 / varargin_1;
  if (!(varargin_1 >= 0.0)) {
    b.set_size(1, 0);
  } else {
    delta1 = std::floor(varargin_1);
    b.set_size(1, static_cast<int>(delta1));
    if (static_cast<int>(delta1) >= 1) {
      nx = static_cast<int>(delta1) - 1;
      b[static_cast<int>(delta1) - 1] = d2_tmp;
      if (b.size(1) >= 2) {
        b[0] = 0.0;
        if (b.size(1) >= 3) {
          if ((-d2_tmp == 0.0) && (static_cast<int>(delta1) > 2)) {
            d2scaled =
                d2_tmp / (static_cast<double>(static_cast<int>(delta1)) - 1.0);
            for (k = 2; k <= nx; k++) {
              b[k - 1] = static_cast<double>(
                             ((k << 1) - static_cast<int>(delta1)) - 1) *
                         d2scaled;
            }
            if ((static_cast<int>(delta1) & 1) == 1) {
              b[static_cast<int>(delta1) >> 1] = 0.0;
            }
          } else if ((d2_tmp < 0.0) &&
                     (std::abs(d2_tmp) > 8.9884656743115785E+307)) {
            delta1 = d2_tmp / (static_cast<double>(b.size(1)) - 1.0);
            i = b.size(1);
            for (k = 0; k <= i - 3; k++) {
              b[k + 1] = delta1 * (static_cast<double>(k) + 1.0);
            }
          } else {
            delta1 = d2_tmp / (static_cast<double>(b.size(1)) - 1.0);
            i = b.size(1);
            for (k = 0; k <= i - 3; k++) {
              b[k + 1] = (static_cast<double>(k) + 1.0) * delta1;
            }
          }
        }
      }
    }
  }
  if (std::isnan(varargin_1) || std::isinf(varargin_1)) {
    delta1 = rtNaN;
  } else if (varargin_1 == 0.0) {
    delta1 = 0.0;
  } else {
    delta1 = std::fmod(varargin_1, 2.0);
    if (delta1 == 0.0) {
      delta1 = 0.0;
    } else if (varargin_1 < 0.0) {
      delta1 += 2.0;
    }
  }
  if (delta1 == 1.0) {
    delta1 = varargin_2 / (2.0 * varargin_1);
    d2scaled = (varargin_1 + 1.0) / 2.0;
    b[static_cast<int>(d2scaled) - 1] = varargin_2 / 2.0 - delta1;
    b[static_cast<int>(d2scaled + 1.0) - 1] = varargin_2 / 2.0 + delta1;
  } else {
    b[static_cast<int>(varargin_1 / 2.0 + 1.0) - 1] = varargin_2 / 2.0;
  }
  b[static_cast<int>(varargin_1) - 1] = d2_tmp;
  wf.set_size(b.size(1));
  k = b.size(1);
  for (i = 0; i < k; i++) {
    wf[i] = b[i];
  }
  if (a.size(1) == 1) {
    firfreqz(1.0 / a[0], varargin_1, h);
  } else {
    nx = static_cast<int>(std::fmax(0.0, 1.0 - static_cast<double>(a.size(1))));
    b_a.set_size(1, a.size(1) + nx);
    k = a.size(1);
    for (i = 0; i < k; i++) {
      b_a[i] = a[i];
    }
    for (i = 0; i < nx; i++) {
      b_a[a.size(1)] = 0.0;
    }
    nx = static_cast<int>(std::fmax(0.0, static_cast<double>(a.size(1)) - 1.0));
    b.set_size(1, nx + 1);
    b[0] = 1.0;
    for (i = 0; i < nx; i++) {
      b[i + 1] = 0.0;
    }
    if (varargin_1 < b_a.size(1)) {
      datawrap(b, varargin_1, a);
      fft(a, num);
      datawrap(b_a, varargin_1, a);
      fft(a, den);
      if (num.size(1) == den.size(1)) {
        h.set_size(num.size(1));
        k = num.size(1);
        for (i = 0; i < k; i++) {
          double ai;
          double ar;
          double bi;
          ar = num[i].re;
          ai = num[i].im;
          d2_tmp = den[i].re;
          bi = den[i].im;
          if (bi == 0.0) {
            if (ai == 0.0) {
              h[i].re = ar / d2_tmp;
              h[i].im = 0.0;
            } else if (ar == 0.0) {
              h[i].re = 0.0;
              h[i].im = ai / d2_tmp;
            } else {
              h[i].re = ar / d2_tmp;
              h[i].im = ai / d2_tmp;
            }
          } else if (d2_tmp == 0.0) {
            if (ar == 0.0) {
              h[i].re = ai / bi;
              h[i].im = 0.0;
            } else if (ai == 0.0) {
              h[i].re = 0.0;
              h[i].im = -(ar / bi);
            } else {
              h[i].re = ai / bi;
              h[i].im = -(ar / bi);
            }
          } else {
            double brm;
            brm = std::abs(d2_tmp);
            delta1 = std::abs(bi);
            if (brm > delta1) {
              d2scaled = bi / d2_tmp;
              delta1 = d2_tmp + d2scaled * bi;
              h[i].re = (ar + d2scaled * ai) / delta1;
              h[i].im = (ai - d2scaled * ar) / delta1;
            } else if (delta1 == brm) {
              if (d2_tmp > 0.0) {
                d2scaled = 0.5;
              } else {
                d2scaled = -0.5;
              }
              if (bi > 0.0) {
                delta1 = 0.5;
              } else {
                delta1 = -0.5;
              }
              h[i].re = (ar * d2scaled + ai * delta1) / brm;
              h[i].im = (ai * d2scaled - ar * delta1) / brm;
            } else {
              d2scaled = d2_tmp / bi;
              delta1 = bi + d2scaled * d2_tmp;
              h[i].re = (d2scaled * ar + ai) / delta1;
              h[i].im = (d2scaled * ai - ar) / delta1;
            }
          }
        }
      } else {
        binary_expand_op(h, num, den);
      }
    } else {
      fft(b, varargin_1, num);
      fft(b_a, varargin_1, den);
      if (num.size(1) == den.size(1)) {
        h.set_size(num.size(1));
        k = num.size(1);
        for (i = 0; i < k; i++) {
          double ai;
          double ar;
          double bi;
          ar = num[i].re;
          ai = num[i].im;
          d2_tmp = den[i].re;
          bi = den[i].im;
          if (bi == 0.0) {
            if (ai == 0.0) {
              h[i].re = ar / d2_tmp;
              h[i].im = 0.0;
            } else if (ar == 0.0) {
              h[i].re = 0.0;
              h[i].im = ai / d2_tmp;
            } else {
              h[i].re = ar / d2_tmp;
              h[i].im = ai / d2_tmp;
            }
          } else if (d2_tmp == 0.0) {
            if (ar == 0.0) {
              h[i].re = ai / bi;
              h[i].im = 0.0;
            } else if (ai == 0.0) {
              h[i].re = 0.0;
              h[i].im = -(ar / bi);
            } else {
              h[i].re = ai / bi;
              h[i].im = -(ar / bi);
            }
          } else {
            double brm;
            brm = std::abs(d2_tmp);
            delta1 = std::abs(bi);
            if (brm > delta1) {
              d2scaled = bi / d2_tmp;
              delta1 = d2_tmp + d2scaled * bi;
              h[i].re = (ar + d2scaled * ai) / delta1;
              h[i].im = (ai - d2scaled * ar) / delta1;
            } else if (delta1 == brm) {
              if (d2_tmp > 0.0) {
                d2scaled = 0.5;
              } else {
                d2scaled = -0.5;
              }
              if (bi > 0.0) {
                delta1 = 0.5;
              } else {
                delta1 = -0.5;
              }
              h[i].re = (ar * d2scaled + ai * delta1) / brm;
              h[i].im = (ai * d2scaled - ar * delta1) / brm;
            } else {
              d2scaled = d2_tmp / bi;
              delta1 = bi + d2scaled * d2_tmp;
              h[i].re = (d2scaled * ar + ai) / delta1;
              h[i].im = (d2scaled * ai - ar) / delta1;
            }
          }
        }
      } else {
        binary_expand_op(h, num, den);
      }
    }
  }
  if (varargin_1 < h.size(0)) {
    if (varargin_1 < 1.0) {
      nx = 0;
    } else {
      nx = static_cast<int>(varargin_1);
    }
    h.set_size(nx);
  }
  nx = h.size(0);
  Sxx_unscaled.set_size(h.size(0));
  for (k = 0; k < nx; k++) {
    Sxx_unscaled[k] = rt_hypotd_snf(h[k].re, h[k].im);
  }
  k = Sxx_unscaled.size(0);
  for (i = 0; i < k; i++) {
    delta1 = Sxx_unscaled[i];
    Sxx_unscaled[i] = delta1 * delta1;
  }
  a__1.set_size(Sxx_unscaled.size(0));
  if (Sxx_unscaled.size(0) != 0) {
    nx = (Sxx_unscaled.size(0) != 1);
    i = Sxx_unscaled.size(0) - 1;
    for (k = 0; k <= i; k++) {
      a__1[k] = v * Sxx_unscaled[nx * k];
    }
  }
  if (std::isnan(varargin_1) || std::isinf(varargin_1)) {
    delta1 = rtNaN;
  } else if (varargin_1 == 0.0) {
    delta1 = 0.0;
  } else {
    delta1 = std::fmod(varargin_1, 2.0);
    if (delta1 == 0.0) {
      delta1 = 0.0;
    } else if (varargin_1 < 0.0) {
      delta1 += 2.0;
    }
  }
  if (delta1 == 1.0) {
    delta1 = (varargin_1 + 1.0) / 2.0;
    if (std::isnan(delta1)) {
      b.set_size(1, 1);
      b[0] = rtNaN;
    } else if (delta1 < 1.0) {
      b.set_size(1, 0);
    } else {
      b.set_size(1, static_cast<int>(delta1 - 1.0) + 1);
      k = static_cast<int>(delta1 - 1.0);
      for (i = 0; i <= k; i++) {
        b[i] = static_cast<double>(i) + 1.0;
      }
    }
    Sxx_unscaled.set_size(b.size(1));
    k = b.size(1);
    for (i = 0; i < k; i++) {
      Sxx_unscaled[i] = a__1[static_cast<int>(b[i]) - 1];
    }
    if (b.size(1) < 2) {
      i = 0;
      nx = 0;
    } else {
      i = 1;
      nx = b.size(1);
    }
    k = nx - i;
    varargout_1.set_size(k + 1);
    varargout_1[0] = a__1[static_cast<int>(b[0]) - 1];
    for (nx = 0; nx < k; nx++) {
      varargout_1[nx + 1] = 2.0 * Sxx_unscaled[i + nx];
    }
  } else {
    delta1 = varargin_1 / 2.0 + 1.0;
    if (std::isnan(delta1)) {
      b.set_size(1, 1);
      b[0] = rtNaN;
    } else if (delta1 < 1.0) {
      b.set_size(1, 0);
    } else {
      b.set_size(1, static_cast<int>(delta1 - 1.0) + 1);
      k = static_cast<int>(delta1 - 1.0);
      for (i = 0; i <= k; i++) {
        b[i] = static_cast<double>(i) + 1.0;
      }
    }
    Sxx_unscaled.set_size(b.size(1));
    k = b.size(1);
    for (i = 0; i < k; i++) {
      Sxx_unscaled[i] = a__1[static_cast<int>(b[i]) - 1];
    }
    if (b.size(1) - 1 < 2) {
      i = 0;
      nx = 1;
    } else {
      i = 1;
      nx = b.size(1);
    }
    k = nx - i;
    varargout_1.set_size(k + 1);
    varargout_1[0] = a__1[static_cast<int>(b[0]) - 1];
    for (nx = 0; nx <= k - 2; nx++) {
      varargout_1[nx + 1] = 2.0 * Sxx_unscaled[i + nx];
    }
    varargout_1[k] = a__1[static_cast<int>(b[b.size(1) - 1]) - 1];
  }
  if (!std::isnan(varargin_2)) {
    k = varargout_1.size(0);
    for (i = 0; i < k; i++) {
      varargout_1[i] = varargout_1[i] / varargin_2;
    }
  } else {
    k = varargout_1.size(0);
    for (i = 0; i < k; i++) {
      varargout_1[i] = varargout_1[i] / 6.2831853071795862;
    }
  }
  varargout_2.set_size(b.size(1));
  k = b.size(1);
  for (i = 0; i < k; i++) {
    varargout_2[i] = wf[static_cast<int>(b[i]) - 1];
  }
}

} // namespace coder

//
// File trailer for pburg.cpp
//
// [EOF]
//
